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We present the simplest nuclear energy density functional (NEDF) to date, determined by only 4 
significant phenomenological parameters, yet capable of fitting measured nuclear masses with better 
accuracy than the Bethe-Weizsacker mass formula, while also describing density structures (charge 
radii, neutron skins etc.) and time-dependent phenomena (induced fission, giant resonances, low 
energy nuclear collisions, etc.). The 4 significant parameters are necessary to describe bulk nuclear 
properties (binding energies and charge radii); an additional 2 to 3 parameters have little influence on 
the bulk nuclear properties, but allow independent control of the density dependence of the symmetry 
energy, excitation energy of isovector excitations, and the Thomas-Reiche-Kuhn sum rule. This 
Hohenberg-Kohn-style of density functional theory (DFT) successfully realizes Weizsacker’s ideas 
and provides a computationally tractable model for a variety of static nuclear properties and dynamics, 
from finite nuclei to neutron stars, where it will also provide a new insight into the physics of the 
r-process, nucleosynthesis, neutron star crust structure, and neutron star mergers. This new NEDF 
clearly separates the bulk geometric properties — volume, surface, symmetry, and Coulomb energies 
which amount to ~ 8MeV per nucleon or up to ~ 2000 MeV per nucleus for heavy nuclei - from finer 
details related to shell effects, pairing, isospin breaking, etc. which contribute at most a few MeV for 
the entire nucleus. Thus it provides a systematic framework for organizing various contributions to the 
NEDF. Measured and calculated physical observables - i.e. symmetry and saturation properties, the 
neutron matter equation of state, and the frequency of giant dipole resonances - lead directly to new 
terms, not considered in current NEDF parameterizations. 
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I. INTRODUCTION 

Calculating nuclear masses, nuclear matter and charge dis¬ 
tributions, and dynamics in nuclear systems remains one of 
the most challenging problems in quantum many-body theory. 
Almost a century ago, Aston ( [1920 1 realized that a nucleus 
is not quite the sum of its parts, leading Eddington ( |1920 1 to 
correctly conjecture a link between nuclear masses, the con¬ 
version of hydrogen into heavier elements, and the energy 
radiated by the stars. When quantum mechanics was first ap¬ 
plied to many-body systems, Weizsacker ( 1935] > proposed that 
an energy density approach could be effective for calculating 
nuclear binding energies. This was the first instance of an 
energy density functional being applied in nuclear physics, 
with the fundamentals of DFT laid several decades later (Drei 


zler and Gross |1990 Hohenberg and Kohn 1 964] |Kohn and 


Sham 1965 i. Bethe and Bacher (1936) further developed 

Weizsiicker's ideas and introduced the nuclear mass formula 
(known as the Bethe-Weizsacker formula) for the ground state 
energies of nuclei with A = N + Z nucleons (N neutrons and Z 
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15.46 

16.71 

22.84 

0 

0.698 

0 

0 

3.30 

15.48 

16.76 

22.88 

0 

0.699 

0 

12.29 

3.17 

15.31 

17.74 

24.94 
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0.765 

-0.667 

0 

2.64 

15.34 

17.78 

24.99 

-22.33 

0.765 

-0.653 

11.46 

2.50 


Table I Parameters and the energy rms of the mass formulas Eqs. 0 
or (2]>, with or without the contribution Eq. © for 2375 nuclei from 
| Audi et al. |(j2012| and Wang et al.\{20\2). (All quantities express in 
MeV.) 


protons): 



N 


E{N,Z) — a v A + a s A~^ + a C ~+ a i -^—— • (1) 

Unlike electrons in atoms, nuclei are saturating systems with 
a nearly constant interior density, yielding the terms above: a 
volume energy, a surface tension, a non-extensive Coulomb 
energy, and a symmetry energy that favors similar numbers 
of protons and neutrons. As shown in the first row of Ta¬ 
ble H] these four terms alone fit the latest evaluated nuclear 
masses ( Audi et al. 2012[ Wang et al. 20121 with a root- 
mean-square (rms) error of Xe « 3 MeV per nucleus, where 
Xe = Y,\En,z ~ E(N,Z)\ 2 /Ne and we fit the Ne = 2375 mea¬ 
sured (not extrapolated) nuclear masses of nuclei with A > 16 
from |Audi et al.\\2Q\2) and [Wang et al.\ ( |20 1 2) i. This is a 
remarkable result: the nuclear binding energy of heavy nuclei 
can reach 2000 MeV, hence the errors are at the sub-percent 
level. 

A slightly better fit is obtained using a mass formula with 
surface corrections terms to the symmetry and Coulomb ener¬ 
gies, as well as odd-even staggering (due to pairing): 


E{N,Z) — a v A +a s A / 3 +ac ^ + a' c ^/t, 


Figure 1 (Color online) The differences £ e xp — Eth in MeVs between 
the evaluated ground state energies (Audi et al. 2012, Wang et al. 
|2012[ l of 2375 nuclei with A >16 and fitted with the 6-parameter mass 
formula Eq. 0 and A = 0. One can easily identify the emergence of 
closed shells for protons and neutrons. 



300 


(. N-Z ) 2 ,((V-Z) 2 

+ Cl1 A +Q1 A 4 /3 +ZL 

(2a) 

Figure 2 (Color online) The binding energy per nucleon B/A = 
\E(N,Z)\/A and the Coulomb, surface and symmetry energy per 
nucleon in Eq. |2|l for the measured 2375 nuclei with A > 16 (Audi 



| et al. 2012} Wang etal., 20121. 

| — (5A~*/ 2 even-even nuclei, 

A = < 0 odd nuclei, 

[ (5A~*/ 2 odd-odd nuclei. 

(2b) 

the geometry and of the periodic trajectories in cavities. As 
early as 1911, Weyl( 1911 1912a|b|c| 1913, 1915 1950) and 
others related the wave eigenstate density in boxes of various 


This pairing contribution is significantly smaller than the oth¬ 
ers, with an amplitude « 12MeV/A 1 / 2 ; it is also smaller than 
contributions arising from shell-correction energies, changing 
the rms error Xe by about at most 100 keV to 150 keV. The 
results of this fit are shown in Table U with the overall residuals 
displayed in Fig. IT] The magnitudes of the various terms are 
compared in Fig. 0 While the volume, surface, and Coulomb 
contributions are clearly dominant to the nuclear binding en¬ 
ergy, the symmetry energy contributions are roughly at the 
level of 10% at most. 

In parallel, a number of properties of many-fermion systems 
were understood mathematically by tying together the roles of 


shapes and boundary conditions to the geometrical shape of 
the box ( |Baltes and Hilf| |1976[ |Brack and BhadurT] |1997[ 


Kac 1966 


Waechter 


1972). In a manner very similar to the 
nuclear mass formula Eq. 0, such formulas can be applied 
to saturating systems, relating the ground state energy to the 
volume (V), surface area (A), and mean curvature radius R of 
the many-particle system: 


E = ayV + a$S + a^R - 


(3) 


The similarity to the nuclear mass formula Eq. 0 becomes 
apparent after relating the volume to the particle number p = 
A/V ss const. (The Coulomb repulsion energy should to be 





































































3 


evaluated separately in a straightforward manner, due to the 
long-range character of the interaction.) The ground state 
energy can thus be rewritten in terms of particle number A 
(here for only one kind of particles) 

E = b v A + b s A 2 / 3 + b R A 1/3 + .... (4) 


The coefficient by is the energy per particle in infinite matter 
and as represents the surface tension. These types of expansion 
are essentially classical in character, with Planck’s constant 
playing no explicit role, and their accuracy for many-fermion 
system is limited by the lack of quantum effects (referred often 
as shell effects). It appears that for nuclei the mass formula 
Eq. |2]) is about as good as one can achieve without introducing 
the quantum effects. 

|Balian and Bloch| ( |1970||197l||l972fr observed that quantum 
states in a finite system can be quite accurately reproduced 
by quantizing the periodic classical trajectories. Combining 
the idea of geometrical quantization with the Thomas-Fermi 
model, the Pauli principle, and copious empirical evidence that 
strongly interacting fermionic systems share many similari- 


ties with non-interacting systems ( 

Haxel et al. , 1949, Landau , 

1956 1957; Mayer, 1949 1950a 

b; Migdal 1967), one can 


quite accurately construct the single-particle density of states 
and binding energies as a function of the particle number, even¬ 
tually correcting this by the shape of the system. 


The single-particle density of states p (e) in a given potential 
has a smooth and an oscillating components: 

p(e) = Ptf(£) + Posc(£), (5a) 

Posc(e) = £flpo(e)sin -Hfcoj (5b) 


where the sum is performed over classical periodic orbits (PO) 
(diameter, triangles, squares, etc.), apo(e) is the stability ampli¬ 
tude, Spo(e) the action, and (j )po the Maslov index of each orbit 
at the energy elBalian and Bloch 1970 1971 1| 1 972] Brack and 


Bhaduri 1997 Nishioka et al. 1 990[ > . The single-particle den 
sity of states in the Thomas-Fermi approximation Ptf (Baltes 


land Hilf||1976l|Brack and Bhaduri|[TW7l[Kacl[T966l|Waechter[ 

dependence on the size and shape of the system, and leads to 
Eqs. ([3| and ([4]) for a square-well potential. At the same time, 
the nature of the periodic orbits also depends on the size and 
shape of the single-particle potential. 

Knowing p(e), one can calculate the particle number A 
and shell-corrections (SC) energy E$c = E — Ej p for a many- 
fermion system by integrating up to the chemical potential 
P: 


ft 1 rv 

A=l p(e)de , Esc = / ep osc (e)de. (6) 

The theory of periodic orbits and structure of these shell- 
corrections has been studied extensively. For example, in a 
3-dimensional spherical cavity, quantum effects can be repro¬ 
duced by including only triangular and square orbits (|Balian| 


land Bloch| [19701 [TWT} [19721 |Brack and Bhaduri| 1 1997] |NIshl 


ioka et al. 1990} >. The emergence of magic numbers, and 
the role of the shapes of many-fermion systems have been 
tested in theory and validated against experimental results 
in fermion systems with up to 3000 electrons (|Brack 1993 


de Heer 1993 Pedersen et al. 19911. In particular, in atomic 


clusters the emergence of the super-shells has been predicted 
theoretically ( |Brack|[T993)|Bulgac and Lewenkopf|[l993||Nish] 
ioka et al. 1990) and confirmed experimentally (|de Heer [1993} 


Pedersen et al. |1991]>, though nuclei are too small to exhibit 


the emergence of super-shells. 

In nuclear physics a similar line of inquiry is encapsu¬ 
lated in the method of shell-corrections, developed by Struti- 


nsky (Strutinsky, 1966 1967 1968) and others (Bohr and 

Mottelson 1998], Brack et al. , 1972 1985 ; Moller et al. 

1995; 

Moller et al. 

2012, Myers and Swiatecki 1974 1990 

1991 

1996 1969 

1966 Ring and Schuck 2004 Strutinsky and 

Magner 197 

5), which shows that p(e) has a well defined 


dependence on the particle number. The smooth part of the 
density of states is well described by the Thomas-Fermi ap¬ 
proximation (and alternatively by the smoothing procedure 
introduced by Strutinsky), the leading terms of which are the 
volume (~ A), surface (~ A 2 / 3 ), Coulomb (~ Z 2 /A 1 / 3 ), and 
symmetry energy (~ (N — Z) 2 /A) contributions encoded in the 
Bethe-Weizsacker mass formula 0 - The oscillating part is 
dominated by the nuclear shape and the shell effects from the 
periodic orbits, with an amplitude dependent on the particle 


number as A 1 / 6 (Strutinsky and Magner 


1976). 


This separation of p(e) into the smooth and oscillating 
parts (|5a|) is a general characteristic of the many fermion 


and Mottelson 1998; Brack et al. 1972, 1985; Moller et al. 

1995 

Moller et al. 

2012 

Myers and Swiatecki 

1974 

1990 

1991 

1996 

1969 

1966. 

ling and Sc 

tuck 5004 

Strutinsky 

1966 

1967 

s 

O' 

00 

Strutinsky and IV 

agner 1976) and self- 


2003; 


Delaroche et al. 


consistent approaches (Bender et al. 

2010} Goriely and Capote 2014| Gorie y et al. 2009 2013 Ko- 

rtelainen et al. | 201 4b|> lead to the same conclusions about the 


various contributions described above, and agree with experi¬ 
mental data (Lunney et al\ [2003] ). This suggests the question: 
To what order can one expand the density of states in powers 
of the particle numbers and periodic orbits? 

There is a reasonable consensus that, beyond the leading con¬ 
tributions from the periodic orbits and shell-corrections, such 
an expansion fails due to the effects of quantum chaos - i.e. 
contributions from classically chaotic trajectories through the 
many-body phase space ( jBohigas et al.\ |T993| ). Stable periodic 
orbits provide the strongest shell effects in quantum systems, 
evident for example in the magic numbers (see e.g. Fig. [4}. 
Unstable periodic orbits also produce shell effects, but typi¬ 
cally with smaller weights. In contrast, chaotic orbits appear 
to produce irregular oscillations in the single-particle density 
of states with a rather small amplitude. Various estimates sug- 


per nucleus ( 

Aberg 

2002 

Barea et al. 

2005 

Bohigas and 

Leboeuf 2002a b; Hirsch et al. 2005 2004; Molinari and Wei- 
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demniiller 2006 2004 Olofsson et al. 2006| 2008] ), noticeably 
smaller than shell effects contributions due to periodic orbits 
and deformations, which are of the order of several MeVs. 

The effect of periodic orbits is not limited to finite systems: 
the Casimir energy in quantum field theory (|Casimir 1948 


Klimchitskaya et al. J2009), critical phenomena (Fisher anc 


de Gennes 1978 Hank e. et al. 1998), and strongly interacting 


infinite inhomogeneous systems, e.g. nuclear pasta phase in 


neutron stars (Bulgac and Magierski 

2001 

2002 Bulgac et al. 

2005 

2006 Magierski et al. 

2003 

Magierski and Bulgac 

2004b 

Magierski and Heenen 

2002b 

|Yu et al. 2000), can also 


be explained and calculated to high precision by evaluating the 
contributions from periodic orbits. This method has become 
the standard for evaluating the Casimir energy in a variety of 


fields i 

Bordag and Pirozhenko 

2010 Canaguier-Durand et al. 

2010|Kjraham 

2014; Rahi et al. 

2009 Schaden 

2010 

). 


It is somewhat surprising that shell effects from periodic 
orbits appear at the same level as deformation effects in the 
energy of nuclear systems. Naively one might expect the defor¬ 
mation energy to be controlled by the surface area of a saturat¬ 
ing system, and thus to contribute as a correction to the surface 
term in nuclear mass formulas like Eqs. 0 and ([2}. However, 
the deformation energy in nuclei has a quantum nature, and 
is determined by a delicate interplay between energy changes 
induced by the changes in surface area and the shell effects. 
A similar behavior has been observed in the case of atomic 
clusters with up to 3000 electrons (Bulgac and Lewenkopf 


1993| . leading to the leveling of the peaks, which one would 


otherwise expect in the absence of deformation), leaving in 
place only the large negative shell-corrections for the magic 
spherical systems, as seen in Fig. [4]in case of nuclei. 

The shape stability of a many-fermion system is controlled 
by the single-particle level density at the Fermi level. In an 
open-shell system this level density is high; the system can 
thus deform quite easily and single-particle levels can rear¬ 
range until the level density is low enough to render the system 
stable. The stabilization process of the nuclear deformation 
in the ground state is analogous to the Jahn-Teller effect in 
polyatomic molecules (Jahn and Teller 1937), where the high 


Our goal is to generate a phenomenological nuclear energy 
density functional (NEDF) with the minimal number of phys¬ 
ically motivated parameters required to describe static bulk 
nuclear properties. This new NEDF may also describe nuclear 
dynamics in real time. 


II. STATIC PROPERTIES 

The lesson from our brief historical review is that, since nuclei 
are saturating systems with a rather well defined saturation 
density, the bulk of the nuclear binding energy should be fixed 
by the geometry of the nuclei (volume, surface area, curvature 
radius) to sub-percent accuracy. As demonstrated in Table [T] 
the accuracy of the mass formulas Eqs. <[T]) and 0 - which 
both lack shell effects, deformation, spin-orbit effects, pairing, 
etc. - suggests that such a nuclear energy density functional 
(NEDF) should be capable of describing at a similar level of 
accuracy both the nuclear binding energies, and the proton 
and neutron matter density distribution. Therefore, we might 
reasonably expect that a NEDF will also describe the nuclear 
charge radii, for which there is a large amount of accumulated 
data (Angeli and Marinova, 2013) >. Quantum effects enter at 
the level of a few MeVs per nucleus, reducing the rms energy 


error Xe from around 3 MeV to about 0.5 MeV (Moller et al. 


|1995||Mbller et a/.~||2012[ ), and are most pronounced for magic 
or semi-magic nuclei, see Fig. [I] The shell corrections will be 
neglected in this round, but can be added through the Strutinsky 
procedure I Bohr and Motte|sonJJ1998| Brack et al. 1972 1985 


Moller et al.\\\ 995||Moller et al. \ |20 1 2] |Myers and Swiatecki 


1974,1990 

1991 

1996 

1969 

1966; Ring and Schuck 

2004 

Strutinsky 

1966 

1967 

1968; 

Strutinsky and Magner 

(976). 


degeneracy of the ground state is lifted by the deformation of 
the system. This mechanism leads to new “magic numbers” 
in deformed systems as Strutinsky discussed in his seminal 
papers (Strutinsky 1966 ; 1967 1968]). The increase in surface 


A number of corrective terms might be considered to im¬ 
prove the accuracy of the nuclear mass formulas Eqs. 0 
and 0. For example, in the Coulomb term, one might re¬ 
place Z 2 with Z(Z — 1) to correctly count the number of proton 
pairs, and one might add an additional term proportional to Z 
to account for the Coulomb exchange interaction and screen¬ 
ing (Bulgac and Shaginyan 1996) 1. Motivated by Eq. ([4]) 
micrhf also consider including terms nronorfional to 4 1 / 3 


one 
3 and 


area and the energy penalty incurred (deformation energy) is 
canceled to a large extent by the shell-corrections (due to pe¬ 
riodic orbits in the deformed potential), unless the system is 
“magic” or “semi-magic”. This cancellation between deforma¬ 
tion energy and shell effects suggests that open-shell systems 
should be easier to deform than magic systems. This property 
is consistent with the character of the residuals remaining after 
the fit of the nuclear binding energies with Bethe-Weizsacker 
formulas like Eqs. (JTJ) and ([2]) and shown in Fig. 0 (see also 
Fig-0- The largest residuals appear as large (negative) spikes 
at the shell closures for spherical nuclei with magic numbers 
of either protons or/and neutrons, while the expected (positive) 
peaks in between magic numbers are flattened. 


might also consider including terms proportional to A 
A 0 . The symmetry energy terms might also be “corrected” by 
replacing (N — Z) 2 /4 with T(T + 1) where T = \N — Z|/2. Fi¬ 
nally, one might introduce an additional correction to account 
for the Wigner energy \N — Z|, which appears as a cusp in 
the nuclear binding energies as a function of N — Z (basically 
only for nuclei with |N — Z\ <2). 

However, including these corrections lead to very small im¬ 
provements in the energy rms Xe beyond the value 2.64 MeV 
obtained with the main terms of Eq. 0- These corrections 
are eclipsed by the shell effects as seen in Fig. 0 From the 
nature of the residuals E exp — E t h in Fig. 0- sharp negative 
spikes at the magic numbers, but roughly constant fluctuations 
in between - one can conclude that mass formulas of the type 
Eq. (0 do encode the role of the nuclear deformation. Shell 
effects are likely responsible for most of these residuals, and 
we shall discuss how to include these toward the end of this 
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paper. A sufficiently accurate theory of nuclear masses may 
even aim to include contributions arising from quantum chaos. 

We will describe a NEDF, which depends on the smallest 
number of phenomenological parameters to account for all the 
contributions in the nuclear mass formulas Eqs. ([TJ and ([ 2 ), 
with the exception of the even-odd staggering due to pairing 
effects. First we relate these parameters to various physical 
quantities relevant for nuclear physics. For a large nucleus, the 
Coulomb energy can be used to estimate the saturation density 
Po by approximating the nucleus as a uniformly charged sphere 
with Ec = 3Z 2 e 2 /5R = acZ 2 /A 1 / 3 , where R = roA 1 / 3 and ro ~ 
1 , 2 fm is a nuclear length scale: 

3 3e 2 

Po = - — 3 , where r 0 = -—. (7a) 

T’/CT'q 1(2 

One can further estimate the ground-state energy of infinite 
nuclear matter per nucleon Eq, the nuclear surface tension a, 
and their dependence on the isospin (N — Z)/2: 


E(N,Z) 
£ ° = A 


= a v + a 1 


(. N-Z ) 2 
A 2 ’ 


a 


— Cl s T Clj 


(N-Z) 2 

A 2 


(7b) 

(7c) 


Finally, one can relate the value of the coefficient a' c (or of the 
alternative coefficient of the contribution a' E Z 2 /A to the mass 
formula (Myers and Swiatecki 19691) with the nuclear surface 
diffuseness. 

One thus expects a NEDF as accurate as the mass formula 
to contain not more than 5 or 6 significant parameters. As 
we shall see, such a functional does exist, and performs com¬ 
parably well with as few as 4 parameters however. That a 
functional dependent on such a small number of phenomeno¬ 
logical parameters can go far beyond the capabilities of mass 
formulae to describe density distributions and dynamics is tmly 
remarkable. 


A. Form of the NEDF 

We postulate a NEDF with three main contributions, 
which significantly improves on the Weizsacker’s original 
idea (|Weizsacker[|1935[): 


4[P//■ Pp —^kin[P«jPp] T $c\pnTpp\ "t“ ^tnt[Pn;Pp]- (8) 

As we shall now discuss, the first two terms are well moti¬ 
vated by a semi-classical expansion and electrodynamics and 
they have no free parameters. All of the phenomenological 
parameters occur in the interaction piece alone. 


Kinetic Terms: The kinetic energy density follows from a semi 
classical expansion of the non-interacting system ([Brack and 


Bhaduri 1997 Dreizler and Gross 1990), and is expressed in 


terms of the neutron/proton number densities p np and masses 


^kin [PniPp\ — 

h 2 
, 2 m t 


E 

z=n.p " 


Vp 


1/2 


' (37T 2 ) 2 / 3 Pt /3 


(9) 


Fligher order corrections to this are considered in Eq. ( f]~3j ) 
below. The factor of 1/9 can be derived for smoothly varying 
densities, and should be compared with the factor of unity 
originally suggested by Weizsacker (1935 1 . 


Coulomb Terms: The direct Coulomb energy and exchange 
contribution in the Slater approximation are: 

Sc(P,uP p ] = \v c (r)p ch (r) - , (10a) 

Vc(r)=|dV^j, (10b) 

where e is the proton charge, and p c h is the charge density, 
which is obtained from the proton and neutron densities by 
convolving with the appropriate charge form factors: 

Pch = G E * p n +G P e * Pp. (10c) 


The charge form factors are determined experimentally, and 
we approximate the Fourier transforms of the form fac¬ 
tors with the dipol e term for the proton , G P E ( Q ) « (1 + 
e 2 /0.71GeV 2 )- 2 |Perdrisat et al.\ 2007), and G n F (Q) « 
a( 1 + Q 2 r\ / 12)- 2 - a( 1 + Q Z G /]2)^hh r 2 = r 2 vg ± 
( r 2 )/2a , (r 2 ) = —0.1147(35) fm 2 , = 0.856(32) fm, and 

a = 0.115(20) (Gentile and Crawford, 2011 1 . Including these 
form factors does not significantly improve the mass fits, but 
improves somewhat the fit of the charge radii. In principle, 
one might allow the coefficient of the Coulomb exchange term 
to vary; this is done, for example, in atomic physics to obtain 
better estimates of the Coulomb exchange energy. We find, 
however, that fitting the nuclear binding energies leads with 
high accuracy to the same coefficient presented in Eq. ( | 1 0a| >, 
so we leave it fixed and do not include this as a parameter in 
our model. 

We require our energy density functional to be an isoscalar 
and include no isospin breaking terms other than those due 
to the neutron-proton mass difference and the Coulomb in¬ 
teraction. Additional isospin violation due to up and down 


quark mass differences and electromagnetic effects (Miller 


et al. | "1990 2006) beyond these two contributions are much 


smaller and are partly responsible mainly for the Nolen- 
Schiffer anomaly (Nolen and Sc hitter |1969) , to which the 
screening of the Coulomb exchange also contributes at a com¬ 
parable level (Bulgac and Shaginyan; i 1996] ). 
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Interaction Terms: All the undetermined parameters of the 
model appear in the interaction terms. We parameterize these 
as 


t &mt(Pn,Pp') — (V o) V. ~ 

T=„,p 


Vp 


1/2 


( 11 a) 


+ ^entrain (Pni Pp) j 
7=0 


g‘ j (p)=a j p 5/3 + b j p 2 + c j p 1/3 , (lib) 

where p is the total density, and (5 is the asymmetry: 


P — Pn + Pp> 


Pn Pp 

Pn T Pp 


(lie) 


We include a correction to the gradient term in ^in so that the 
overall coefficient of the gradient term is rj , and up to 9 pa¬ 
rameters aj, bj , and Cj for j £ { 0 , 1 , 2 } describing the equation 
of state for homogeneous nuclear matter. Three of these (for 
j = 2 ) will be fixed by the equation of state of neutron matter 
determined in ab initio calculations (see section pTCj i. Three of 
the remaining six parameters (a o, c\, and a combination of a \ 
and b\) are found to be only marginally significant at the level 
of changing the energy rms by §Xe < 0.1 MeV, so that in the 
end we shall be left with only 4 significant parameters: rj, bo, 
Co, and a combination of a\ and b\. 

We shall later include an entrainment term Entrain with a 
single parameter a, the form and meaning of which will be 
elucidated in section IV see Eq. ( |40} . We find this term to 
have very little significance for the static properties of nuclei, 
but we will show that it has a pronounced effect on isovector 
dynamics, such as the giant dipole resonance (GDR) mode. 
For the moment while we discuss static properties we set it to 
zero a = 0. One more parameter, |, will be discussed when 


we will discuss the computation of shell effects in section IV.B 


The equations that determine the equilibrium densities of 
a nucleus are obtained by minimizing the energy of a given 
nucleus E(N,Z) = f d 3 rS’\p n ,p p ] with respect to the densities, 
while constraining the total numbers of neutrons N and protons 
Z with two chemical potentials 


n _2 1/2 . TT 1/2 1/2 

~ri- -V p r ' +U x p t = PrPr , 

2m z 

dS[P n ,Pp\ r 

Ur = - = -— , for T e {, n,p}. 

o Pr 


( 12 a) 

(12b) 


B. Gradient Terms 


Weizsacker (|1935) i originally introduced the term proportional 
to gradients of the densities |Vp^^ | 2 with a value of /] = 1 
that was later shown to be valid only if the density has small 


amplitude rapid oscillations ( 

Brack and Bhaduri 

1997 Drei- 

zler and Gross , 1990; Jones and Gunnarsson 1 

989 

). It was 


later rigorously proven that a semi-classical expansion of the 


non-interacting fermions in the extended Thomas-Fermi ap¬ 
proximation (the limit of a slowly varying external potential) 
yields the value 7 ]tf = 1 /9, which defines the lowest order 
gradient contributions in rip ln ([9]) ( Brackand^hadun^ 1997 


|Dreizler and Grossl|1990| [Jones and Gunnarsson 1989| l. We 

treat the coefficient rj as a phenomenological parameter since 
gradient terms can also be generated by interactions (|Gebre- 


|mariam et «/.||20f0]|Negele and Vautherin||1972||1975 1. Fit¬ 

ting the nuclear masses yields values of T) close to 0.5, roughly 
half-way between the semi-classical and Weizsacker values. 

A semi-classical expansion of non-interacting 
fermions ( |Brack and Bhaduri] 1 1 997[ [Dreizler and Gross|p~990| l 
yields the following higher order corrections to that we 
have not included in Eq. ®: 


^(Pn,Pp) ~ ^ 2wt 810 (3^2)2/3^)> 


/(P) = P 


= oV3 


Vp 

p 


27 (YR) 2 YlR + 3 ( v2 P 


(13) 

2 ' 


This type of correction has been studied in nuclear physics 
and shown to lead to quite accurate estimates of the kinetic 
energy density within the extended Thomas-Fermi approxi¬ 
mation ( |Brack et~ai\ |1985||1976[|Brack and Bhaduri} \\991) . 
As within a DFT, such terms can also arise due to the finite 
range of the interactions in a matter similar to some Skyrme 
interactions (Gebremariam et al. 2010[ Negele and Vautherin 


1972 1975 i. However, these terms - even with adjustable 

parameters - do not significantly change the quality of the 
mass fits, so we do not consider them in our main analysis. 
Including them perturbatively in the fit, however, does improve 
the fit of the charge radii. For example, fitting the overall 
coefficient reduce the charge radii residual Xr (see details in 
Section III I from Xr ~ 0.14fm to Xr ~ 0.09fm. Fitting each 
of the three terms independently further reduces the residuals 
to Xr ~ 0.06 fm. We do not include such fourth-order terms 
in our functional as they can lead to a complex behavior of 
the emerging equation for the densities, which can be difficult 
to rationalize. (See, for example, the analysis of fourth order 
differential equations arising in case of non-local potentials 
by |Bulgac| ( |1988| >.) Higher order gradient corrections than 
Eq. ( [T3 | i lead to an unphysical behavior of the densities in the 
classically forbidden regions. The semi-classical expansion 
has an asymptotic character (Jones and Gunnarsson |1989) . It 
is known that corrections beyond second order do not always 
lead to improvements (Jones and Gunn arsson);1989| . 

We point out one more property that will be discussed in 


more detail in section IV a value of 77 = 1/4 corresponds to a 
dynamical theory of superfluid neutron and proton pairs. One 
might naively have thought of T] « 0.5 as an effective nucleon 
pair mass /w e ff ~ 2 m, but this leaves the potential U T wrong 
by a factor of 2. The parameter r\ may simply be thought of 
as a way to control the falloff of the densities in the surface 
region where the interaction effects are still strong. One should 
not expect to obtain the correct asymptotic behavior of the 
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densities far from the nucleus, where the interactions become 
vanishingly small. The main reason is that DFT is by definition 
an integral approach, which minimizes globally the energy and 
not particular local quantities. 

A gradient term alone of the form |Vp| 2 instead of 
(Vp */ 2 | 2 = | Vp j 2 /4p leads to unphysical density profiles, with 
a discontinuity in Vp at a finite radius, beyond which the den¬ 
sity vanishes exactly. 



p,! /3 [fm '] 


Figure 3 (Color online) The quantum Monte Carlo (QMC) re¬ 
sults of Wlazlowski et al. (2014) for the interaction energy per 
neutron displayed as the ratio E n /E F Q defined in Eq. < fT3b> (with 
/3 = 1), where e FC = 3h 2 (3 k 2 p n ) 2 / 2 / ( 10m„), If a„ = 0 in Eq. ( | 15b| ), 
the ratio £j nt /Efg would tend to 0 for p„ —► 0. For densities 
1 /3 

Pn |°aw| < 1 (where = 18.9fm is the .s-wave neuttron-neutron 
scattering length) the leading order correction to the kinetic energy 
density per particle contribution would be instead linear in density 
4 7th Cl]yNp n /l7l n . 


C. Infinite Nuclear and Neutron Matter 


In infinite homogeneous nuclear matter, as might be found in 
a neutron star for example, the gradient and Coulomb terms 
vanish (charge neutrality is maintained by a background of elec¬ 
trons). Neglecting the small neutron-proton mass difference 
m n ss m p = m, the functional acquires a simpler form: 


, 3/i 2 (37T 2 ) 2 / 3 5/3 5 / 3 , 

^(PniPp) ~ \Qm n "F Pp ) 

2 


+ E ( a jP 5/3 + b jP 2 + c jp 1/3 ) P 2J > ( 14 ) 
f=o 


This portion of the functional is essentially an expansion in 
powers of the Fermi momenta k F : k n p = ( 3 n 2 p ntP ) 1 ^ with 
three terms only k Fl k p and k 7 F , an ubiquitous power expansion 
in many-body perturbation theory. It is also supported by fitting 


the neutron matter equation of state (p p = 0, [5 = 1): 

£n{pn) = ^—(37t 2 p n ) 2/3 +£m{p n ), (15a) 

t'mt ( Pn ) — &nPrJ ~\~b n p n -\- C n pn ■ (15b) 

where 

a n =ao+a\+a 2 = — 32.6MeVfm 2 , (15c) 

b n =b 0 +b 1 +b 2 =-n5A Me Vfm 3 , (15d) 

c„=co + ci+c 2 = 109.1 MeVfm 4 . (15e) 

These values have been obtained by fitting the neutron matter 
equation of state as calculated with QMC including up to N 3 LO 
two-body and up to N 2 LO three-body interactions from chiral 
perturbation theory ((Wlazlowski et al. 1|2014), as shown in 
Fig.0 Thus, we can include directly the fixed j = 2 parameters 
ci2, b2, and C2 from the values of a n , h n , and c„ describing the 
QMC results without adding additional free parameters to the 


III 


adding the quartic /j 4 


NEDF. 

As we shall discuss in section 
(j = 2) terms does not significantly affect the quality of the 
mass fits since most nuclei realize small asymmetries ft < 0.25. 
This demonstrates an important point: nuclear masses do not 
constrain the quartic terms /l 4 . Higher powers thus provide 
a direct (and independent) handle on the equation of state of 
neutron matter. 

At this time we do not have an equally accurate QMC cal¬ 
culation of nuclear matter with varying isospin composition, 
so we must rely instead on a phenomenological approach. Our 
main assumption is that we can describe both the isoscalar 
(j = 0) and isovector (j = 1, /j 2 ) parts of the nuclear equation 
of state in a similar fashion to the equation of state of pure neu¬ 
tron matter, using the same three powers of the corresponding 
Fermi momenta. This approach differs from typical Skyrme- 
like parameterizations, which include terms with higher powers 
of densities (e.g. p X/ 4 arising from tp type of terms, where T 
is kinetic energy density). 

The terms ajp 5 / 3 are somewhat unexpected and they are not 
included in Skyrme-like parameterizations. Tondeur (1978) 
introduced only a term a i (without theoretical justification), 
but it makes sense to include the other a ; - for several reasons. 
The QMC calculations of Gandolfi et al. ( 2015| l; |Gezerlis and] 
|Carlson| ( |2010| l; and | Wlazlowski et al. (2014), see Fig. [3] are 
consistent with the existence of a non-vanishing parameter 
a n in the neutron equation of state, which implies that a„ = 
L/=o a i / 0. They also appear naturally in the case of the 
unitary Fermi gas (Zwerger 2012) , a fact confirmed to high 
precision in many experiments. The unitary Fermi gas is a 
system of two species of fermions, interacting with an ,v-wave 
interaction with zero range and infinite scattering length. In 
response to the Many-Body X challenge posed by Bertsch in 
1999, Baker ( 1999) showed that the system was stable. The 
energy density of the unitary Fermi gas scales exactly like the 
kinetic energy density of a free Fermi gas S p 5 / 3 . Since both 
neutron and protons have similar ,v-wave interaction properties, 
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one expects the nuclear energy density to behave somewhat 
like the unitary Fermi gas. 

Although the energy density of the unitary Fermi gas scales 
as the kinetic energy, this is not necessarily due to a mass renor¬ 
malization as one might naively suspect and the answer is not 
a clear cut one unless additional information is obtained. The 


QMC calculations of the single quasi-particle dispersion {Carl¬ 


son and Reddy 2005) and the calculation of the spectral weight 
function ( Magierski et al.\ |2011||2009| l, both find almost the 
same effective mass « m in the unitary Fermi gas, very close 
to unity. Flowever, this does not preclude the interpretation 
of some part of this term as arising from the kinetic energy 
density T as is the case in the unitary Femi gas too ( jBulgacj 
|2007|[20l3)|Bulgac et a/.||2012a| i. The QMC calculations are 
not of enough accuracy to either completely exclude or confirm 
an effective mass different from unity. 

Furthermore, QMC calculations of dilute neutron matter 
at both zero and finite temperatures point to a very large 
pairing gap A k, 0.25 Ef (Gezerlis and Carlson 2010) and 
A = 2.87; (jAbelmcrSeki 2009a|b, Wlazlowski and Magier¬ 


ski 2010 Wlazlowski and Magierski 2011), where £p and 


T c are the Fermi energy and critical temperature for the super¬ 
fluid to normal transition. In comparison, Bardeen-Cooper- 
Schrieffer (BCS) superfluids or mean-field approaches lead 
Awl J6T C , indicating that neutron matter is not a simple BCS 
superfluid (Magierski et al. 2011|[2009). This is further sup¬ 


ported by evidence for a pseudogap at temperatures above T c 
where Cooper pairs and a significant depletion of the density 
of states still exist even in the absence of long-range order. The 
properties of neutron matter are thus closely related to those 
of high T c superconductors. The common approach of adding 
pairing correlations within an Hartree-Fock-Bogoliubov (HFB) 
or BCS approximation must thus be considered only qualita¬ 
tively valid, and properly tuned DFT approach is required for 


2012a: 2011 

2012b Bulgac and Yoon. 2009; Bulgac and Yu 

2002 2003, 

tetcu et al. ||2014 20111 Wlazlowski et al .[ 2015 


Yu and Bulgac[|2003a|b[ ). 


III. FITTING MASSES AND CHARGE RADII 

We fit our NEDFs to the Ne = 2375 measured (not extrapo¬ 
lated) nuclear masses with A > 16 from ( [Audi et ~al\ |2012[ 
Wang et al.\ 2012). We also consider the N r = 883 match¬ 
ing charge nuclear radii from ( |Angeli and Marinova| |2013[ ) 
with Xr = L|<5r| 2 /7V,- When we include the charge radii in 
the fit, we minimize the following quantity Xe/{ 3MeV) 2 + 
Xr / (0.05 fm) 2 which roughly equalizes the weight of the mass 
and radii contributions in the fit. 

At this point, we have 7 parameters in our NEDF: 77 , ao,i, 
7>o,i, and co.i (the j = 2 parameters are fixed by the neutron 
matter equation of state). In addition, we include by hand the 
conventional even-odd staggering Eq. ( |2b| ) with a coefficient 8 
to describe pairing correlations, even though this has very little 
significance in the fits. We consider the following fits: 
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Figure 4 (Color online) The blue dots show the results obtained using 
fit NEDF-1 with Xe = 2.58MeV, while the red dots are the results 
of the fits using nuclear mass formula Eq. 0 - with Xe = 2.52MeV. 
When compared against each other, the rms energy deviation between 
the two fits is Xe = 0.36MeV. Thus, NEDF-1 essentially reproduces 
the nuclear mass formula Eq. 0 - The lower panel is a rendering of 
calculated ground state energy differences between the experimental 
measured values and NEDF-1 results, in which one can see clearly 
the magic numbers separately for neutrons and protons. 


NEDF-0: A six parameter least-squares fit of the Ne = 2375 
nuclear masses (Audi et al. 2012| Wang et al. 2012 1 
including rj, bo, co, at, b\, and 8 but setting the nucleon 
charge form factors Eq. ( |10c[ > G E = 1 and G n E = 0. 


NEDF-1: The same as NEDF-0, but including the measured 
charge form factors. Comparing with NEDF-0 we see 
that the electric form factors are not significant for the 
overall mass fits, but slightly impact the charge radii at 
the 0.01 fm level (for the reduced x? )- 
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NEDF 

11 

ao 

bo 

co 

a\ 

b\ 

Cl 

a 2 

b 2 

C2 

8 

Xe 

[MeV] 

Xr 

[fm] 

0 

0.4719 

0 

-3.15166 

2.12378 

1.048 

-0.610 

0 

0 

0 

0 

0.3246 

2.59 

0.145 

1 

0.4742 

0 

-3.13778 

2.10995 

0.981 

-0.544 

0 

0 

0 

0 

0.3250 

2.58 

0.135 

2 

0.4743 

0 

-3.14595 

2.11873 

0.961 

-0.521 

0 

0 

0 

0 

0 

2.71 

0.140 

lr 

0.4807 

0 

-2.98351 

1.94501 

1.087 

-0.668 

0 

0 

0 

0 

0.3330 

2.71 

0.051 

3 

0.4800 

-0.088 

-2.95408 

2.01459 

1.000 

-0.549 

-0.017 

0 

0 

0 

0.3234 

2.58 

0.139 

3n 

0.4739 

-0.061 

-3.00477 

2.03797 

1.565 

-1.367 

0.181 

-1.765 

3.8811 

-1.9731 

0.3279 

2.57 

0.133 

3 nr 

0.4815 

-0.061 

-2.86708 

1.89073 

1.563 

-1.347 

0.169 

-1.763 

3.7239 

-1.8138 

0.3529 

2.67 

0.050 

E 

0.4885 

0 

-3.14903 

2.11957 

0.277 

0.277 

0 

0 

0 

0 

0.3177 

2.64 

0.129 

Er 

0.4957 

0 

-3.00642 

1.96664 

0.264 

0.264 

0 

0 

0 

0 

0.3600 

2.74 

0.051 

En 

0.4866 

0 

-3.15157 

2.12271 

0.272 

0.272 

0 

-0.533 

2.3889 

-1.8763 

0.3192 

2.62 

0.133 

Enr 

0.4970 

0 

-3.00488 

1.96492 

0.260 

0.260 

0 

-0.521 

2.2540 

-1.7185 

0.3545 

2.74 

0.051 


Table II Dimensionless fit parameters and residuals for the various NEDFs. Parameters have been scaled by appropriate powers of py = 
0.15fm -3 and e F = £ (37T 2 p 0 /2) 2 /3 = 35.2941961307MeV: Sj = ajp^ 3 /e F ,bj = bjp 0 /e F , and cj = cjp^/ep. 


NEDF-2: The same as NEDF-1, but without the pairing pa¬ 
rameter 5 = 0. Comparing with NEDF-1 we see that 
odd-event staggering is also relatively insignificant at 
the level of 0.1 MeV per nucleus. This is consistent with 
the results from the mass formulas in Table Q] 

NEDF-lr: The same as NEDF-1, but including the N r = 883 
charge radii into the fit. We see that there is significant 
room to improve the description of the charge radii with¬ 
out significantly degrading the mass fits. The charge 
radii residuals for NEDF-1 and NEDF-lr are plotted in 

Kg- ID 

NEDF-3: The same as NEDF-1, but with all 8 parameters, in¬ 
cluding no and ci that we omitted from the previous fits. 
In conjunction with the principal component analysis 
shown in Fig. El this fit demonstrates that the terms 
with parameters ao and c\ are insignificant. 

NEDF-3n: The same as NEDF-1, but with all 8 parameters, in¬ 
cluding ao and c i that we omitted from the previous fits, 
and the /l 4 parameters for the terms quartic in isospin, 
constrained by the QMC neutron matter equation of 
state ( [Wlazlowski et «/.~||2014[ ) using Eqs. ( |15b| >. That 
the quality of the fit, isoscalar, and isovector parameters 
change very little, demonstrates that the neutron matter 
equation of state is essentially independent of the nuclear 
masses. 

NEDF-3nr: The same as NEDF-3n but including the charge 
radii as in fit NEDF-lr. That the ao and ci terms are 
insignificant for both masses and radii is also emphasized 
by this fit. 


To test this, we set a\ = b\pl^ where p* = 0.15fm 3 is 

1 /3 

a constant. The combination a\ — b\pj , to which the 
masses are insensitive, allows independent control the 
slope Li of the symmetry energy (see Eq. ( |32e[ >). From 
the fits we see that this same combination also controls 
the neutron skin thicknesses. 

NEDF-Er: The same as NEDF-E but including the charge 
radii as in fit NEDF-lr. 

NEDF-En: This is our main fit. It is the same as NEDF-E 
but includes the /3 4 parameters adjusted to reproduce the 
neutron matter equation of state as in fit NEDF-3n. 

NEDF-Enr: The same as NEDF-En but including the charge 
radii as in fit NEDF-lr. 

These fits are summarized in Table HU with the saturation and 
symmetry properties in Table |Hl| The residuals for fit NEDF-1 
are shown in Fig. [4] and compared with a fit to the nuclear with 
mass formula Eq. Q. 

The reduced Xe for these fits is comparable to that obtained 
using the nuclear mass formulas Eq. (JTJ) with 4 parameters (plus 
5) and Eq. ([2]i with 5 parameters (plus 5). This is consistent 
with our hypothesis that a NEDF for masses should contain no 
more than 5 significant parameters. Note, however, that unlike 
the mass formulas, the NEDF also gives a good description 
of charge radii - for which the mass formula says nothing - 
and provides access to nuclear dynamics as we shall discuss in 
section El 

A. Discussion 


Finally, we consider two parameter sets that represent the 
minimal functionals. Each has 4 significant parameters (and 5, 
which is insignificant): 

NEDF-E: Following the principal component analysis of 

NEDF-3n (discussed below) we find the combination 

1 /3 

a\—b\pj to be only weakly constrained by the mass fit. 


The accuracy of the ground state binding energies obtained 
using this NEDF, see Fig.[4]and Table |II| compares extremely 
well with much more sophisticated self-consistent approaches 
that naturally account for the shell-corrections. 

The UNEDF2 nuclear energy functional introduced by |Ko-j 
rtelainen et al. (2014a|b i has a residual of Xe = 1.95 MeV 
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per nucleus for 555 even-even nuclei. Since the UNEDF2 
functional is formulated in terms of single-nucleon orbitals, 
one would expect it to account for shell-corrections, but this 
functional still displays large discrepancies of the binding en¬ 
ergies of magic nuclei (close to 8 MeV in case of 20 X Pb). The 
UNEDF2 functional depends on 11 parameters, with several 
additional parameters to describe pairing correlations. The 
UNEDF2 functional also leads to a larger pairing coupling 
constant for protons than for neutrons, thus violating isospin 
invariance and the expectation that Coulomb effects reduce 
proton pairing, an expectation which appears to agree with 


experiments (Lesinski et al. 2009 1 . A natural solution might 
be to include in the NEDF a pairing contribution of the form 


^pairing — goiPtnPp) [|V n | + |Vp| ] 

+ 8l(pn,Pp) [N 2 -|Vp| 2 ] \p„-pp], 

gO.\(Pn,Pp) = gO.l (PpiP;i)) 


(16a) 

(16b) 


where v„ p are the anomalous neutron and proton densities 
respectively. Since in measured nuclei one has predominantly 
N>Z, see Fig. [7] a phenomenological analysis that leads to an 
apparent coupling for protons larger than the one for neutrons 
can be reconciled with coupling constants go(PmPp) < 0 and 
8l(j>n,P p) >0 - 

|Baldo etTd. (2013 2008) 1 introduced an energy density func¬ 
tional based on information extracted from QMC calculations 
of neutron and symmetric nuclear matter, and a few additional 
parameters to describe pairing and spin-orbit interaction and 
finite range effects; they assume that no quartic terms in isospin 
/j 4 are present in the NEDF. Depending on how they fit var¬ 
ious subsets of about 580 nuclei, they find Xe ranging from 
1.3 MeV to 2.4 MeV. 


Goriely et al. (2009 20131 have produced over the years a 


series of NEDF in which they obtain an even lower Xe through¬ 
out the entire mass table as low as 0.5 MeV. However this was 
obtained by adding a number of phenomenological corrections, 
a procedure which so far has not been adopted by the other 
practitioners using microscopic approaches. 

Within relativistic mean-field theory (RMFT) of nuclei in 
one of the best parametrizations of the NEDF for nuclear 


masses (Agbemava et al. 2014 Niksic et al. 20111 one 
achieves a Xe between 2 and 3 MeV for even nuclei using 
the AME2012 data set ( |Audi et a/.||20T2l|Wang et a/.|[20l2) i. 


B. Neutron Matter 

The inclusion of the j = 2 terms quartic in /j 4 in fits NEDF-3n 
and NEDF-3rn demonstrates an important point: the equation 
of state of neutron matter has very little impact on the form of 
the NEDF at the level of quadratic isovector contributions (see 
also the previous section). We have attempted to perform a fit 
of the nuclear binding energies by including the quadratic terms 
in /j only, defining the coefficients of the functional through 
the neutron equation of state ( a\ = a n — a o, b\ = b n — bo, and 
ci = c n — co), and allowing 77 , ao, bo, and co to vary. The 
energy rms was at best Xe = 4.36 MeV. 


In measured nuclei, the ratio )3 = (p„ — p p )/p ~ (N — Z)/A 
is \P 1 < 0.25 (with a very small number of exceptions), hence 
nuclear masses are essentially insensitive to the presence of the 
/3 4 terms, as |/31 4 < 1 /256. Only 46 nuclei have \N — Z\/A > 
1 /4 (mostly with A < 50) and only 68 nuclei have N < Z in 
the set we consider. To assess the magnitude of these effects, 
we have evaluated the /i 4 contributions to the nuclear binding 
energies perturbatively, see Fig. [5] This contribution is quite 
small and can be easily overlooked when discussing known 
nuclei, but is crucial in order to correctly reproduce the energy 
of neutron matter. 


Thus, using properties of the neutron matter to constrain 
the form of the NEDF and/or arguing against the inclusion of 


higher powers of (p„ — p p ) ( 

Baldo et al. 2013 

2008 Brown 

2013; Brown and Schwenk 

2014;|Erler et al. 

2013;|Fantina] 

et al. 2014, 

Fayans( 1998; Goriely et al. \ 2013 

Reinhard and 

Nazarewicz 

2010 ) is an ill-advised procedure, and the applica- 


tions of functionals constructed in this manner, in particular to 
star environments, should be regarded with suspicion. 
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Figure 5 (Color online) The contribution to the ground state energies 
of the terms quartic in isospin density 8 E 14 = J rf 3 ^(p )/3 4 , evaluated 
perturbatively with NEDF-1, see Table [Tl] In the lower panel we 
display the ratio (N — Z)/A for the nuclei we have considered. Among 
the 2375 nuclei we have considered, there are 33 nuclei with N = Z, 
78 nuclei with Z > N, and 70 nuclei with \N — Z\/A> 1 /4. 
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NEDF 

Po 

[fm- 3 ] 

-E0 

K 

S 

L 

l 2 

Neutron skin 
208 pb 48 Ca 

[fm] [fm] 

0 

0.136 

15.24 

222.5 

26.8 

34.1 

32.8 

0.082 

0.118 

1 

0.136 

15.22 

222.4 

26.7 

35.9 

34.7 

0.087 

0.123 

2 

0.136 

15.21 

222.2 

26.7 

36.8 

35.6 

0.089 

0.125 

lr 

0.148 

15.48 

227.7 

27.1 

30.9 

29.6 

0.078 

0.116 

3 

0.136 

15.21 

216.5 

26.7 

34.7 

33.4 

0.088 

0.124 

3n 

0.137 

15.20 

218.2 

30.0 

29.3 

16.7 

0.068 

0.107 

3 nr 

0.147 

15.44 

222.9 

31.0 

31.2 

15.5 

0.068 

0.107 

E 

0.136 

15.28 

223.1 

29.7 

68.2 

66.9 

0.159 

0.174 

Er 

0.147 

15.53 

228.1 

30.6 

70.2 

68.9 

0.161 

0.176 

En 

0.136 

15.27 

222.9 

30.1 

29.1 

66.1 

0.152 

0.172 

Enr 

0.147 

15.53 

228.2 

31.1 

31.1 

68.3 

0.156 

0.174 


Table III Saturation, symmetry, and neutron skin properties for the 
various NEDFs. All values in MeV unless otherwise specified. 



C. Saturation and Symmetry Properties 

The isoscalar parameters j = 0 and quadratic isovector param¬ 
eters y = l C/3 2 ) may be directly related to the saturation and 
symmetry properties respectively by expanding the energy per 
nucleon of homogeneous nuclear matter Eq. ( [T4| about the 
symmetric saturation point p„ = p P = po / 2 : 

S[pn ' pp) = £o (p) + £ 2 (p)j3 2 + £4(p)j3 4 + ^(j3 6 ). (17) 

r 

The saturation density po, energy per nucleon £q, and incom¬ 
pressibility Kq are then defined by the minimum £q(Po) = 0, 
and depend only on the j = 0 isoscalar parameters «q, bo, and 
co- Expanding about po in 5 = (p — po)/3po and in powers of 
13 = (p„ — p p )/p, one can define various “local” contributions 
to the symmetry energy .S'i. 4 , its density dependent slope Li. 4 , 
etc.: 


£ o(p)^£o + + C?(5 3 ), 

e 2 (p)=S2-L2d + ^K 2 d 2 + ff{8 3 ), (18) 

£ 4 (p) = S 4 -U8 + ±K 4 8 2 + 0{8 3 ) 


Since we include also quartic terms /l 4 , we must differentiate 
between these local symmetry parameters S 2 , L 2 , etc. and the 
full symmetry parameters defined as the difference between 
symmetric matter and pure neutron matter (see also the discus¬ 
sion of |Lattimerjp014| ): 


S = 


S(po,Q)-S{po/2,po/2) 


L=3p 


dp 


Po 

<f(p, 0 ) 


= 3p 0 £'(po). 


Po 


(19) 

( 20 ) 


where £„p is the neutron equation of state ( | 15 a[ > . Since the 
saturation density po minimizes the energy of symmetric mat¬ 
ter, the slope of the full symmetry energy L at p 0 depends 
only on the equation of state of pure neutron matter. Thus, the 
QMC neutron equation of state alone fixes the global density 


Figure 6 (Color online) Equation of state of symmetric nuclear matter 
(all functionals) and neutron matter (for NEDF-3n, 3nr, En, Enr only) 
in the lower panel and pure neutron matter (upper panel) for NEDF-0, 
1. lr, 2, 3, E, and Er (which do not constrain the neutron equation 
of state and have only |3 2 contributions). The black dots shows the 
QMC results (Wlazto wski el q/.1 |2Dl4). The symmetry energy 5 is 
indicated for the functionals NEDF-Er and Enr. The slope L ~ 30MeV 
is fixed by the neutron matter equation of state alone (if used as 
a constraint, see Eq. (|20|>). In this case the slope L 2 /3po may be 
tuned without significantly affecting the mass fit by a djusting the 


1/3 

-b\pj , see Section 


III.G 


Functionals 


insensitive combination a \ ■ 
with only quadratic isospin contributions (/3 2 ) appear to cross near 
p„ ~ 0.1 fm -3 (upper panel). 


dependence of the symmetry energy L = 3po£/,(po) ~ 30MeV. 
(Small variations in L seen in Table |IH| for functionals fitting 
neutron matter are due to the variations in po from fit to fit.) 

When only /3 2 isospin contributions are included in the func¬ 
tional, our fits to the nuclear binding energies display a fea¬ 
ture reported in other NEDFs: the energy per neutron in pure 
neutron matter appears to be well constrained at a density of 
p„ « 0.1 tin 4 where all functionals cross, see Fig. [ 6 ] However, 
the value for the energy per neutron ss 9 MeV at this point in 
our fits is significantly smaller than the value k, 12.35 MeV 
obtained in QMC calculations of |WlazIowski el al.\ (|2014| . 
This feature is not present when the /3 4 terms are included 
(NEDF-3n, NEDF-3nr, NEDF-En and NEDF-Enr) and the 
QMC results are reproduced without significantly affecting 
the mass fits. The statement often made in the literature (see 
Horowitz et «/.~j( |2014a] ) and references therein) that the value 
of the symmetry energy at p ~ 0.1 fm 1 is well constrained 
by nuclear masses must only be applied to the local expansion 
S 2 at this density, but not to the symmetry energy difference S 
between symmetric and pure neutron matter. 
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Figure 7 (Color online) The proton and neutron driplines evaluated with NEDF-1, NEDF-3n, NEDF-E and NEDF-En compared with those 
predicted by fully self-consistent calculations using SkM*, SkP, SLy4, SV-min, UNEDFO, and UNEDF1 i jErler el q7~||20 1 2} . The 2375 nuclear 
masses from ( Audi et al. 2012, Wang et al.. 20121 are displayed as dots and nuclei unstable with respect to the emission of a single proton or 
neutron are denoted with open circles. NEDF-E leads to a neutron dripline shifted quite a bit away from the other predictions. This is due 
to changing L 2 to a value of the neutron skin thickness in 20X Pb consistent with experiment, see Table 
equation of state and the values for S and L, which appear to control the position of the dripline. 


Ill 


but inconsistent with the neutron 


D. Symmetry Energy and Neutron Skin Thickness 


As shown in Table III the binding energy of nuclear matter and 
the symmetry energy predicted by all NEDFs fits agrees well 
with the value obtained with the mass formula (|2j. Our fits 
generally estimate the slope of the symmetry energy Lo from 
29 MeV to 36 MeV, but NEDF-E, NEDF-En, NEDF-Er, and 
NEDF-Enr demonstrate that this is essentially unconstrained 

by the masses and can be adjusted independently with the 

1/3 

combination a\—b\pj . 

We also compute the neutron skin thickness of 48 Ca and 
208 Pb, for which precision measurements CREX and PREXp] 
are underway, see (Horowitz et al. 2014b 1 for details. For 
NEDF-1 through NEDF-3m, the neutron skin of 208 Pb ranges 
from 0.07 fm to 0.09 fm while the 48 Ca skin ranges from 
0.106 fm to 0.125 fm. The 208 Pb neutron skin appears quite a 


bit thinner than a recent measurement of 0.15(3) fm (Tarbert 


1 Proposals and related information available at http: //hallaweb. j lab 
org/parity/prex 


[eFa7j[[20T4) , but the results of NEDF-E, NEDF-Er, NEDF-En, 

and NEDF-Enr demonstrate that this is also controlled by the 

1/3 

same combination ci\—b\p* as Li, and hence unconstrained 
by the masses. 

Since the slope of the symmetry energy L = Li + L 4 H-« 

30 MeV is fixed by the neutron matter equation of state, requir¬ 
ing a larger value of L 2 ~ 60 MeV to explain the neutron skin 
thickness of 208 Pb also suggests that at least quartic terms are 
required in the functional. 


E. Neutron Drip Line 

It is interesting to compare the limits of nucleon stability pre¬ 
dicted by our NEDFs. In Fig. [7] we compare the proton and 
neutron drip lines obtained with our NEDF-1, 3n, E, and En 
agains the predictions of UNEDFO and UNEDF1, as well as 
those obtained with other Skyrme parametrizations extracted 
from the supplemental data of Erler et al. ( |2012j >. It is a bit 
unexpected that NEDF-1 and NEDF-3n and En have almost 
identical neutron driplines considering that NEDF-1 and E 
have only quadratic isovector contributions, while NEDF-3n 































13 


and En also have quartic terms constrained by the neutron 
matter equation of state. Apparently these quartic terms - and 
therefore the neutron equation of state - play a minor role even 
in nuclei with a very large neutron excess N — Z ss 160 when 
(N — Z)/A ss 0.4 since [(N — Z) /A ] 4 ss 0.0256 is small. 

The proton driplines are in good agreement with other 
NEDFs (up to odd-even effects), even though our NEDFs 
lacks quantum effects. Our NEDFs suggest that the neutron 
dripline may extend somewhat further than for conventional 
functionals. If significantly more nuclei are stable against nu¬ 
cleon decay, this may dramatically impact the astrophysical 
r-process, which is predicted to follow lines of constant separa¬ 
tion energy in close proximity to the neutron dripline |Langake| 
|and Wiescher| |200 1 \ |Meyer| |1989[ > . |Meyer| ( |1989| l considered 
neutron star ejecta as the site of r-process nucleosynthesis, and 
determined that the reaction flow is very close to the dripline. 
Even though his simulations were performed for relatively 
cold matter (recent simulations seem to indicate that the star 
material is somewhat heated (Go riely et «/.| [20TT] |Rosswog| 
| et al. 1 12014| >), it will be interesting to simulate the r-process 
using the present NEDFs. There are almost 9000 nuclei be¬ 
tween the proton and neutron driplines in case of NEDF-1 
and NEDF-3n, most of which might be stable against nucleon 
decay (quantum effects are neglected in this estimate). This 
number is noticeably larger than the estimate 6900 ± 500 ob¬ 
tained in Ref. ( |Erler et ai 1 12012| >. The functionals NEDF-E 
and NEDF-Er lead to a neutron dripline which is shifted much 
further, but when the neutron equation of state is incorporated 
in NEDF-En and NEDF-Enr, the neutron dripline is again very 
close to the previous position given by NEDF-1 and NEDF-3n. 
Unlike the neutron skin thickness, which is controlled by Li, 
the position of the neutron dripline appears to be controlled by 
the full symmetry energy S and its density dependence L. 

The predicted position of the neutron dripline will likely 
affect the structure of the neutron star crust inferred from older 
studies (B aym et al. 1971 1 Bulgac and Magierski| 2001 J2002J 


Magierski and Bulgac 2004a| Magierski et al. J2003 lMagier- 


ski and Heenen 2002a| Magierski and Bulgac |2004b| Negele 


and Vautherin 1973) . The corresponding increase in the neu¬ 

tron skin thickness will also affect the profile and the pinning 
energy of quantized vortices in the neutron star crust ( Avo-| 


gadro et al. 2007 2008 Bulgac et al. | 2013 |Pizzochero et al. 


1997 Pizzochero 2007 2011 Yu and Bulgac||2003c l. 


Fusion cross sections (Adelberger et al. |2011 Gasques 
et al. 1 12005| ) will also be significantly altered, particularly 
in stellar environments where neutron rich nuclei fuse via 
pycnonuclear reactions ( Afanasjev et al. | |2012| |Schram and| 
Koonin| 19901, and where the neutron gas surrounding nuclei 
leads to their swelling ( |Umar etal. 2015| >. A thicker neutron 
skin with further enhance this effect. 


F. Charge Radii 

Using the parameters determined from the mass fits, the NEDF 
also models the neutron and proton densities in the nuclei. 
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Figure 8 (Color online) The difference between the measured <|An-| 
geli and Marinova 2013) and computed charge radii of 883 nuclei 
evaluated with NEDF-1 (blue) and NEDF-lr (red), see Table[I] The 
energy rms corresponding to the blue and red points are 2.58 MeV 
with NEDF-1 and 2.71 MeV with NEDF-lr. The charge radii rms are 
in these cases 0.135fm with NEDF-1 and 0.051 fm with NEDF-lr 
respectively. 



0 2 4 6 8 10 

r [fm] 

Figure 9 (Color online) The calculated proton p p (r) (dotted) and 
charge p c h(r) (solid) densities for 48 Ca (red) and 208 Pb (blue), calcu¬ 
lated with NEDF-En (lower curves) and NEDF-Enr (upper curves) 
compared to charge densities (black) extracted from electron scatter¬ 
ing experiments EE Vries et oT||1987| >. 


allowing us to extract the charge radii for the 883 nuclei in I An- 


geli and Marinova 

2013) that intersect with the 2375 nuclear 

masses we fit from ( 

Audi et al. 2012] Wang et al. 2012 ). With- 


out any further adjustment, the NEDF gives an rms residual 
for these charge radii of % r = 0.135fm. Including them in 
the fit (NEDF-lr, NEDF-3nr, and NEDF-Enr) improves this 
to X> = 0.05 fm without significantly affecting the quality of 
the mass fits. To compare, fully self-consistent RMFT models 
realize Xr from 0.023 fm to 0.041 fm ([Agbemava et al. 2014 


|Niksic et al. | |201 1| ) depending on the set of nuclei examined. 

The differences between the calculated and measured charge 
radii (Angeli and Marinova, 2013| ) are illustrated in Fig. [ 8 ] In 
Fig. M we compare the proton and charge densities of ^Ca 
and ^Pb calculated with NEDF-E and NEDF-Er with the 
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Figure 10 (Color online) The various ellipses show the region in 
the (e 0 ,Po) plane, in which the NEDF parameters can be changed 
and to lead to changes in the residual S%e < 0.2MeV. While the 
equilibrium energy Eq and density po are controlled mainly by the 
combination bo + cp, which is constrained with very high precision, 
the combination bo — cq remains significantly less constraint, see 
Section |HI.G| This aspect allows us to manipulate to a certain degree 
the saturation properties, while affecting the overall fit only slightly. 


the charge densities extracted from electron scattering experi¬ 
ments (De Vries et al. 1987[ l. The calculated charge density 
208 pb bas a slightly larger radius and slightly smaller diffuse¬ 
ness compared to those extracted from data. There may be 
several ways to cure this deficiency. The nuclear saturation 
density is not well constrained by the mass fits, and our values 
appear somewhat low. However, one can constrain the equi¬ 
librium density to a value closer to po = O.lbfm 3 without 
affecting the mass fit greatly, see Fig. [lO] The surface proper¬ 
ties are controlled by the parameter r ), and although a single 
parameter is sufficient to fit the masses, one might introduce a 
density dependence r\ —> rjo + t]iP 1/<3 + • • • or even an isospin 
density dependence to better fit the densities. 

We do not have yet a clear understanding of why the light 
nuclei have slightly smaller radii, while the large nuclei show 
the opposite behavior, and whether this aspect is correlated 
with the nuclear surface diffuseness. Within the present ap¬ 
proach all nuclei are spherical; quantum effects may have a 
noticeable impact on the calculated charge radii, but so far 
we do not know how to disentangle the roles of the bulk and 
shell-correction effects on these values as cleanly as we do in 
the case of binding energies. 

In order to compare our results with the experimental val¬ 
ues of the charge radii, we must correct for the motion of the 
center of mass, but we find these corrections to be small. In 
mean-field-like calculations, the center-of-mass of the system 
moves in the self-consistent potential and one needs to account 
for these fluctuations before comparing with experimental data. 
For small nuclei, it is common to treat these center-of-mass 
fluctuations by assuming that the center-of-mass sits in a har¬ 
monic oscillator potential, with the frequency related to the 


average curvature of self-consistent nuclear mean-field poten¬ 
tial (Lipkin 1958; Negele |1970| l. The calculated charge radius 
would then be smeared by convolving with the gaussian center- 
of-mass wavefunction, and the correction implemented as a 
de-convolution which amounts to the subtraction 


,r?) = jc?rr‘p c lr>-- 2 


3 Ti 


2 Amco 


/ jd : 'rr 2 p c (r)-A 2,/3 1.8fm 2 , 


( 21 ) 

( 22 ) 


where hco ~ 1,85MeV + A _1 / 3 35.5MeV approximates the 
harmonic-oscillator energy (Negelej j 1970 1 . (In the full cal¬ 
culation, the charge form factors are included in a similar 
manner by subtracting the proton and neutron charge radius.) 
A similar correction of jhco should be subtracted from the 
binding energy of the nucleus. This prescription is expected 
to be accurate enough for light nuclei, but it is not clear that 
the harmonic approximation of the mean-field potential will be 
valid for large nuclei, which saturate to approximately po in the 
core. (Due to the effect of the nuclear surface tension the nuclei 
are compressed a bit, but at the same time Coulomb repulsion 
tends to deplete the central charge density.) Here we argue 
that, for large nuclei, this correction due to the center-of-mass 
fluctuation is small enough to be discarded (an argument made 
also by jNegelej(1970| for large nuclei). 

To simplify the following argument, we neglect the proton- 
neutron mass difference; in actual calculations we use the 
experimentally measured nucleon masses. The intrinsic radius 
of a nucleus r mls is the expectation value of the square of the 
distance of the nucleons from the center of mass: 


ns = lE<(^-*) 2 > = lE<^>-<* 2 > 


^ k= 1 


A 

A ^v* 2 } 
A k= 1 

= rl - (R 2 ) 


(23) 


where R = Y*k r k/A is the center-of-mass coordinate and \Jrf n 
is the rms mass radius of the nucleus calculated directly within 
the meanfield theory. This makes it clear that the center-of- 
mass correction - the second term - is a measure of the size of 
the center-of-mass fluctuations 

(* 2 > = X2 (E r ? + E 'V'-/') 

A \k=i k^i / 

The last term in Eq. ( |24| is the value of a two-body observable 
and in Hartree-Fock approximation only the exchange term 
contributes. Using the known expression for the two-body 
density in the Hartree-Fock approximation 

Pi(x,y) = pi(x,x)pi(y,y) -pi(x,y)pi{y,x), (25) 

where p i (x,y) = p i (ri, tfi, X\ , < 72 , Ti) is the one-body den¬ 

sity matrix. We now approximate this using the Thomas-Fermi 
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approximation for the density matrix 


Pi(*))') — 5oi,o 2 &ti, t 2 [ x3 


where kp is the local neutron or proton Fermi momentum and 
s = r\ — V 2 , one can show that 


1 

A 3 






(27) 


Thus, the leading order terms of the direct and exchange contri¬ 
butions cancel exactly leaving ( R 2 ) = which is signif¬ 

icantly smaller than the correction suggested by the harmonic 
oscillator assumption ( Lipkin| 1958, Negele 1970). We there¬ 
fore neglect these correction in our fits as they are insignificant 
in comparison with missing shell effects etc. 

It is instructive to perform a similar calculation for the ki¬ 
netic energy of the center-of-mass 


Erm — 


1 







k= 1 



(28) 


Y.Pk Pi 


(For clarity we have suppressed spin and isospin indices.) Us¬ 
ing again the Thomas-Fermi approximation ( |26| ), and including 
the Weizsacker correction (Brack and Bhaduri 1 997} |Dreizler| 
and Gross 1990| Weizsacker 1935[ ) to the kinetic energy den¬ 
sity in the direct term, we determine that 


2mA 


j d 3 r (Vp 1 / 2 (r) 


2 mrZ 


A ” 1 / 3 


(29) 


where m is a nucleon mass, R = roA 1 / 3 is the nuclear radius, 
and tr/Imr^ ~ 35MeV is of the order of the Fermi energy. 

The arguments presented here on the role of center-of-mass 
fluctuations are based on a Hartree-Fock-like approach, since 
we need to evaluate two-body observables, and strictly speak¬ 
ing within DFT, one does not have access to the two-body 
densities. 


G. Principal Component Analysis 

The parameters listed in Table [II] are highly correlated. To ana¬ 
lyze these, we consider as significant changes 8%e ~ 0.1 MeV 
since this is the typical level of sensitivity of the mass fits. We 
keep the changes relatively small because otherwise the model 
is not well approximated by a quadratic error model if 8%p > 
0.1 MeV. Numerically we find that even 0.1 MeV is too large, 
but yields qualitatively correct information after a full refitting. 
Note that <S(xJ) = (x E + 8 xe) 2 ~Xe = 2 Xe 8xe + (8xe) 2 , 
so we must normalize 8(Xe ) by 2 Xe 0.1 MeV in order to 
consider changes Sxe ~ 0.1 MeV. 




Figure 11 (Color online) Principal component analysis for the NEDF- 
1 fit (top) and the NEDF-3 fit (bottom). Plotted are the components 
of the eigenvectors v„ defining the principal component Eq. pib| > as 
linear combinations of the various dimensionless parameters. From 
this we see that the most-significant component po ~ + c0 which 
fixes the saturation energy to high precision. At the same time the 
component pn ~ (bo — cq in NEDF-1 (and similarly in NEDF-3n) is 
not well constrained. We also see that the least-significant component 
P 5 ~ a\ — b\ is essentially unconstrained. For NEDF-3, we find three 
insensitive components, two of which can be used to set the smallest 
parameters ag = c\ = 0. After removing these, one obtains a similar 
analysis as for NEDF-1 above. 


To compare the parameters in a meaningful way, we must 
make them dimensionless and of order unity. We do this as in 
Table |n]by scaling with appropriate powers of p* = 0.15fm 3 
and Ep = |^(37r 2 p*/2) 2 / 3 = 35.29420MeV, which we take 
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as fixed parameters close to the saturation values: 


a jP 


2/3 


£f 


h J = 


bjP* 

£p 


ajp 


4/3 


£f 


(30) 


(It is important to retain a significant number of digits for 
isoscalar quantities, as it will be come clearer below.) In partic¬ 
ular, we consider the covariance matrix C such that the residual 
deviation is 


——«s r -c- 1 -s = y;&^ 

2Je • 0.1 MeV V A 2 


(31a) 


where 8 is the deviations vector of the dimensionless param¬ 
eters Eq. ( |30| ) from their best fit values as listed in Table [II] 
and we have diagonalized Cv„ = A 2 v„ to obtain the principal 
components p„ 

P„ = v« • (do b 0 ••• rj 8). (31b) 


Since the parameters are of order unity, we may directly con¬ 
sider the A„ as a measure of the errors: changing p„ by A„ will 
affect the fit on the scale of Sxe ~ O.lMeV. Therefore, the 
smaller the value of the parameter A„, the more precisely the 
fit to nuclear masses constrains the value of the corresponding 
linear combination of NEDF parameters. 

The principal components for fits NEDF-1 and NEDF-3 
are shown in Fig. m Both show similar features that can be 
understood in terms of the saturation and symmetry parameters 
Eqs (fl 8 | which we may express as: 


S 

e F 

L 

£f 


— — I +<5o+ b 0 + coj (32a) 

£ F 

0=3 -l-flo-l-jfeo T2co, (32b) 

—pj—= 5 +do - 2c 0 , (32c) 

—— 5 —~ + id.\ + b\ +ci) + (<22 + b 2 + C 2 ); (32d) 

2S 

-b(&i+2ci) +(^2 + 2cb). (32e) 


The most significant component po in both fits is the sum of the 
7 = 0 coefficients do + bo + do which fixes the saturation energy 
£0 Eq. ( |32a[ ) to very high precision. Next are mixtures of r| 
and the symmetry energy S Eq. ( |32d| >, which are correlated 
by the finite size of the nuclei; the latter is the sum of the 
7 = 1 coefficients d\ +b\ +ci (recall that these fits have no p 4 
components). Finally, we have some insignificant parameters, 
including <5. 

In the case of NEDF-3, we see that two parameters are 
completely insignificant. These include do ~ —0.088 and c 1 = 
—0.017. These values are an order of magnitude smaller than 
the other coefficients: hence, the insignificant components can 
be easily removed by setting «o = ci = 0 which we do in most 
of our fits. 


Finally, both plots indicate that a combination of the 7 = 1 
parameters is highly insignificant. Thus, in NEDF-1, the com¬ 
bination b\ — a\ can be given almost any value of order unity 
without changing Xe more than O.lMeV. This is directly 
tested in NEDF-E, NEDF-Er, NEDF-En, and NEDF-Enr where 
we change the sign of b\, setting a\ = b\p* . Indeed, we see 


that Xe changed by about O.lMeV. Notice from Table III 


that the slope of the symmetry energy L 2 changes from about 
30MeV to 70MeV while the other parameters remain about 
the same. This also significantly changes the neutron skin 
thickness, demonstrating a correlation between L 2 and the 
skin thickness, similar to that seen in other mean-field mod¬ 
els (Vinas et al. 2014| >. This is consistent with Eq. ( |32e[ > where 
we see that b\ gives us a direct handle on Li. 


IV. DYNAMICAL PROPERTIES 


The NEDF can be used to study a variety of dynamical phenom¬ 
ena, excitation of nuclei by various external probes and various 
(large amplitude) collective modes, low energy collisions of 
nuclei, induced fission and fusion, etc. A dynamical inter¬ 
pretation of the functional follows by considering a coupled 
hydrodynamic theory of protons and neutrons. 

Consider first a single species with number density p and 
a local momentum potential (j) such that the local fluid mo¬ 
mentum p = mv = V0. The Euler-Lagrange equations for the 
following Lagrangian density 

^=-p(^+^^) 2 ]-W-V^ n ^VP) 2 (33a) 

define the following hydrodynamic equations: 


p+V(pv) = 0 , 

• mv 2 . h 2 V 2 Vp 


=0, 


the latter of which is sometimes expressed as 


mv- 




= 0 . 


(33b) 

(33c) 


(33d) 


Recognizing that the last term in these equations has the form 
of a “quantum pressure” term with a modified Planck’s con¬ 
stant fi = 17 1 / 2 /z, one can check that this hydrodynamic theory 
is equivalent to the following “quantum” theory for a wave- 
function t f/ = v / pexp(i 0 /^) = v / pexp(i 0 ): 

/ Z -\ 

= lihdt + ^-jy-gip), (34a) 

h 2 

ifi\j/= --^-V 2 i// + e?'(p)i//. (34b) 

2m 

Note that with these equations, the action is invariant under the 
Galilean boost to a frame with velocity v: 




fid 


mv 2 t 

mv ■ x — 


2 


(35) 
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Equations ( [34] ) have exactly the same form as a non-linear 
Schrodinger equation that would be used to study the collective 
wavefunction of a superfluid Bose-Einstein condensate (BEC) 
with homogeneous energy density S'(p). The implementation 
of the hydrodynamic equations as a non-linear Schrodinger 
equation allows for stable numerical implementations and may 
be solved efficiently using standard algorithms. Somewhat 
surprisingly this “Schrodinger equation” admits “quantized 
vortices”, albeit with a modified “Planck’s constant” ft. The 
quantization of vortices here is purely a formal artifact since 
it depends on the free parameter I]. Although it might be ap¬ 
pealing to consider Eqs. ( |34| ) as a theory for superfluid dimers 
comprising two protons or two neutrons, the best-fit value of 
77 « 0.5 does not presently admit this interpretation. As men¬ 
tioned earlier, to restore properly quantized vortices requires 
77 = 1/4 and a rescaling of the wavefunction, by normalizing 
it to the total particle number - i.e. 2 p = |i//j 2 - so that the 
factors of 77 =1/4 cancel in Eqs. m yielding a properly 
normalized Schrodinger equation. 

One can similarly recast our complete NEDF in terms of two 
“complex wavefunctions” with a modified “Planck’s constant” 
fi = ^yfjti: 

tyn,p = y/Piwp ex P(iQn,p), m np V np ~ • (36a) 

through the Lagrangian density 

^ = ih{ y/,j y/,, + yl \j/ p ) - S[yn , y p \ (36b) 

where 

fp ff _ 

£[Vn, w P ] = ^- v v4 • V V4 + • v w+ 

4“ ^entrain ( tyn ? typ) 4" -‘h (Pn iPp)i (36c) 


Ignoring for a moment the entrainment term (^entrain ( Yn, Yp), 
dynamics can be implemented using the following coupled 
non-linear “Schrodinger” equations: 


~ dy/,, p 
ill y 


dt 
Un.p 


fl 2 

^ ^ Yn,p 4“ Un : p Yn,p> 

~ u >n.p 

d $h(PiuPp) 

d Pn,p 


(37a) 

(37b) 


A. Entrainment 


As mentioned above, this theory is covariant under Galilean 
boosts to a frame with velocity v (Eq. (|35j)). Since both protons 
and neutrons experience the same boost v np —>• v„ p + v , we 
can introduce an additional Galilean invariant term in the two- 
component system proportional to \v„ — v p \ 2 , obtaining the 
collective kinetic energy: 


m„p n vl 

2 


m pPp v l 

2 


! t/ ^h!^lpPnPp 7 '// 

"i oc ~ 

2mp 



(38) 


The parameter a controls the entrainment of the neutron and 
proton fluids, and the form has been chosen so that the term 
vanishes when either density vanishes. 

As with the quadratic terms above which enter the functional 
as Vi/Vp) 2 , the entrainment term must be constructed from 
the wavefunctions rather than from the velocities to avoid 
singularities when the density vanishes (i.e. in the cores of 
quantized vortices or in the tails of the nucleus). This can be 
done by noting that the following “currents” transform under 
Galilean boosts as 


Jn.p — ' ft V4; tyn p ■ 


J n,p ^ Jn,p 4” HhupPn.pV • (39) 


and the homogeneous energy-density (no gradients) is 


Hence, we define the Galilean invariant entrainment term: 
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E 
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+ E ( a jP 5/3 + a jP 6/3 + a iP 1/3 ) 


4/3 


Pn Pp 

P 


2 j 


(36d) 


The representation of the densities in the form ( |36a| ) is similar 
in spirit to the decomposition of the single-particle density 
matrix championed in the theory of large amplitude collective 


motion p —> e ix pe ix (Baranger and Veneroni 

00 

Brink 

\et a/.l[1976||Ring and Schuck 

2004 

). The energy 

of a nucleus 


is thus also uniquely separated into the collective kinetic energy 
(the terms explicitly depending on Yn,p) and the internal energy 
depends only on the “shape” of a nucleus, with no internal 
quasi-particle excitations. 


^entrain (Y»,Yp) = 

t HhilllpPnPp 
2mpo 
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ITlnPn 


= ce¬ 


ll 1 

2 mp 


1/2 
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1/2 

mn 


nipPp 
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Yp^Yn~ l/2 Yn^Yi 


For static properties, this term has the form 


(40) 


^entrain (P« > Pp) — 


= a- 


2 mp 


1/2 


1/2 


p; ,! vp,! /2 - '2//p,! /2 vpP 


(41) 


and leads to a coupling between the neutron and proton density 
gradients Vp„ - Vp ;) . Bodmer (1960) introduced this type of 


term, and it also appears in various Skyrme parameterizations 
of the NEDF. By varying a from —0.5 to 0.5, Xe changed 
by at most 15 keV. The significant effect of this term is seen 
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however in the dynamics, where the motion of one fluid will 
drag along the other, affecting strongly the excitation energies 
of isovector modes such as the GDRs. 

Entrainment (the Andreev-Bashkin effect) was predicted by 


Andreev and Bashkin (19751 to occur in superfluid mixtures 
of J He and 4 He, and is rather surprising at first sight, since 
superfluids are expected to flow without resistance. In particu¬ 
lar, one might have expected that if somehow one would bring 
into motion only one superfluid component, superfluidity will 
have the consequence that the other component remains at rest. 
The entrainment term ( [40] ) is indeed dissipationless, and thus 
it does not violate superfluidity, but allows the motion of one 
superfluid to influence (entrain) the other. It is natural to expect 
a similar phenomenon to arise in nuclei, where proton and neu¬ 
tron (super)fluids can coexist. Entrainment should also plays a 
role in neutron stars and has been studied intermitently since 


and Haensel 2006 Gusakov and Haensel, 

2005; Vardanyan 

and Sedrakyan 

1981 Volovik et al. 1975 

). The formalism 


describing these systems is called three fluid hydrodynamics - 
two superfluids and one normal component - and is generaliza¬ 
tion of Landau’s two fluid hydrodynamic phenomenological 
model of superfluids at finite temperatures below, but the crit¬ 
ical temperature. Since in nuclear systems both neutron and 
proton subsystems can have a superfluid and a normal compo¬ 
nents at finite temperatures, and since the normal components 
can move independently in isovector modes, a proper general¬ 
ization of the superfluid dynamics to nuclear systems would 
be a four fluid hydrodynamics, with two superfluid and two 
normal components, thus a somewhat more complex system 
than the superfluid mixtures considered so far in literature. 


B. Shell Effects 

Shell-corrections have a typical amplitude of several Me Vs. 
Even though they represent a relatively small correction to the 
total energy of the system, the can have a major effect on low 
energy nuclear dynamics. Their magnitude is controlled by 
the spin-orbit interaction, the pairing effects, and the nuclear 
shape. There are several prescriptions one can use to compute 
shell corrections, but the idea behind them is the same. The 
shell-corrections are a function E sc ({fcy}) of single-particle 
energies £& 


hrYr,k = £x,kYr,k, (42) 

where i f/ T k are the single-particle wave functions for the single¬ 
particle Hamiltonian h x for proton and neutrons T = n,p. This 
Hamiltonian is the functional derivative of the energy density 


(here a = 0) 


^\p„,P P ]= E ^|vw| 2 
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i\ r 


9 ) 2in r 
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1/2 


T=n.p 

' E + $c(pn,Pp) 

j =0 

spin—orbit + pairing contributions. 


(43) 


The main difference between this energy density and Eq. ([8]) 
is that now the Yr,k are the single-particle neutron and proton 
wave functions. 

There is an ambiguity is separating the contribution of the 
gradient terms and other forms of this energy density can be 
considered. For example, one can introduce an additional 
parameter (which might have a density and even an isospin 
dependence) and define the energy density &\p n ,p P \ as: 


^\Pn,Pp)= E (l+S)£-\VY T ,k\ 2 

k y T=n,p T 


E £-!(3" 2 ) 2/3 Pr /3 


x=a, P 2m t 


i\ r 
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T=n,p \ 

- E S ’j(P')P 2j + S 'c(Pn;Pp) 

7=0 

spin—orbit + pairing contributions. 


(44) 


The introduction of such a parameter | requires adding addi¬ 
tional terms to restore the Galilean invariance of this energy 
density functional (Engel et at 1975) (strictly speaking the 
Galilean invariance of the intrinsic energy of a nucleus), and 
its introduction will also lead to a modification of the Thomas- 
Reiche-Kuhn sum rule for GDR. The semiclassical limit of 
both Eqs. ( |43| and ( |44| are identical by construction. 

As in the case of the unitary Femi gas ( |Bulgac||2007||2013] 
Bulgac ~et al.\ 2012a|>, without additional information form 


either the spectral weight function or from spectroscopy, the 
interpretation of the terms a jp 5 ^ as either being the semiclas¬ 
sical limit of the appropriate kinetic energy densities or of 
interaction energy densities is ambiguous. This ambiguity is 
related to the discussion above. 

The equations of motion ( |37a| ) have to be augmented now 
with the contribution to the potential arising from E sc ({as 
follows: 
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Figure 12 (Color online) Induced fission of 238 U excited with an 
initial quadrupole and octupole potential velocity 6„(r) = 0 p (r) = 
<72(2 z 2 —x 2 —y 2 ) + q 3 (5z 2 — 3)z starting from the ground state, see 
Eqs. {36a}. The symmetric case (left) has q^ = 0 and an excitation en¬ 
ergy of 7.68 MeV (initial collective kinetic energy). The asymmetric 
case (right) has an excitation energy of 10.33 MeV (initial collective 
kinetic energy) and the final charges, neutron numbers, and atomic 
masses of the two fragments are Z\ = 49.3, A) = 78.4, A\ = 127.7, 
Z 2 = 42.7, N 2 = 67.6, and A 2 = 110.3. Both fragments emerge with a 
significant amount of collective excitation energy and several excited 
multipolarites. (Note that the interval between the first three frames is 
about 5 times the interval between the lower frames.) 


The dynamics described by Eqs. {45} is in spirit similar to the 


constrained density TDHF approach developed by Oberacker 


|and Umarj{2015} and |Umar and Qberacker|{2006| ). The nuclear 
systems will evolve in time, while only the collective degrees 
of freedom are active, in exactly the same manner as the theory 


of Adiabatic TDHF (ATDHF) was envisioned (Baranger and 


Veneroni] |1978| Brink et al. 1976| Ring and Schuck| 2004} . 

The present approach is formulated directly in terms of (all) rel¬ 
evant collective degrees of freedom and there is no difficulty in 
separating the degrees of freedom into intrinsic and collective. 
The definition of the collective degrees of freedom and their 
proper separation from the intrinsic ones is a problem prac¬ 
tically unsurmountable in the usual theory of LACM {Dang| 
et al ., |2000} . The dynamics described by Eqs. {45} , or their 
simplified form Eqs. ( | 1 2a| ), is by default isentropic. There is no 
need and no difficulty in constructing either the inertia tensor 
or the potential energy surface, as the system naturally evolves 
from one point to another in the collective phase space. 


C. Induced Symmetric and Asymmetric Fission 


In Fig. [12] we show the evolution of a 238 U nucleus starting 


from the ground state with two different initial isoscalar veloc¬ 
ity kicks and no shell-corrections, leading to symmetric and 
asymmetric fission respectively. The fragments emerge with 
a significant excitation energy and deformation, clearly not 
in the corresponding ground states. The excitation energy of 
the nucleus has two components, the collective flow given by 
Eq. {38} and the deformation and Coulomb energies given by 
Eq.{8} respectively. The total energy is naturally conserved. 
The results presented in these figures have been obtained by 
solving the full three-dimensional (3D) time-dependent equa¬ 
tions {34} , on a 64 3 spatial lattice. This takes about 30 minutes 
on a MacBook Pro using a MATLAB program about 250 line 
long (including the generation of the ground state and of the 
movie and related graphics). 


D. Giant Isovector Resonances 

Equations {34} have also been used to extract the behavior 
of the GDRs as a function of the atomic mass A. As a func¬ 
tion of the atomic mass, the NEDFs without entrainment yield 
IKOgdr ~ 65... 70/A 1 / 3 MeV, which is too low (Piotr Magier- 
ski, private communication). Even though the entrainment 
term plays a small role when computing the ground state en¬ 
ergies, it plays an important role in determining the excitation 
energies of isovector modes. 

Let us consider for simplicity a nucleus with N = Z, in which 
we will neglect the Coulomb effects and the neutron-proton 
mass difference m « m n « m p . In this case in the ground 
state p„(r) = p p (r) = po(r)/2 = |<//o(r)| 2 and the density are 
determined form the nonlinear equation 

fi 2 

— V~Wo + &o(Po)Vo = )Wo- 
2m z 


( 46 ) 
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Considering small amplitude isovector modes p niP (r,t) = 
|y/b(r) ± y/\((r,t)\ 2 , where y/o(r) = i//q( r), one can derive the 
equation: 


ihW = -(r7 + a)—V 2 vn+2#o(Po)v^i 

Am 

iQ £i{Po), , *s s 2 /V(vwi)\ 

+8—-—(vri + v^)-a—vr 0 V- —^— . 

Po w V r 0 / 


(47) 


For a > 0 the effective mass is lowered (first term) and the 
stiffness increased (last term), which results in an increased 
value of the ScBgdr- Thus, the entrainment term through the 
parameter a provides a direct handle on isovector dynam¬ 
ics without significantly impacting the static properties of the 
NEDFs. One last remark: Eqs. ( |34| > will only provide a good 
description of the collective modes and their corresponding 
transition strengths. 


V. CONCLUSIONS 


masses and the neutron matter equation of state are essentially 
uncoupled. 

These NEDF are suitable to study complex collective motion 
such nuclear fission, low energy collision of nuclei, their fu¬ 
sion, interactions with various external time-dependent probes, 
structure and dynamics of the neutron star crust, nucleosyn¬ 
thesis, etc. in a computationally efficient manner in order to 
better gain intuition and understanding about nuclear dynam¬ 
ics. Finally, an entrainment term with parameter a controls the 
energies of the isovector nuclear excitation modes and together 
with an additional parameter 8, one can control the value of the 
Thomas-Reiche-Kuhn sum rule. 

This NEDF contains several new elements with respect to 
commonly used Skyrme-like density functional theories: 

• Terms proportional to p 5 / 3 , similar to those found in the 
study of the unitary Fermi gas. 

• The gradient terms in Eq. ([8]) proportional to 
p|V P y 2 | 2 = ? 7 1Vp T | 2 /4p T orto aV p,! /2 • V p p /2 . 


The nuclear energy density functional (NEDF) developed in 
this work contains at most 7 significant parameters, each clearly 
related to specific properties of nuclei. In particular, to globally 
fit the measured ground state energies of 2375 nuclei, function¬ 
als NEDF-En and NEDF-Enr can achieve an rms residual of 
Xe ~ 2.6 MeV to 2.7 MeV with only 4 significant fit parame¬ 
ters (and an insignificant pairing parameter 8). 

Therefore NEDF-En and NEDF-Enr reproduce the behavior 
of the best nuclear mass formula, with fewer free parameters, 
and yet have many advantages; unlike nuclear mass formulas, 
these NEDFs also model the neutron and proton densities in 
the nuclei, allowing one to extract the charge radii. Fitting the 
masses alone gives an rms residual of x? = 0.13 fm for the cor¬ 
responding 883 measured charge radii, while including them 
in the fit improves this to Xr = 0.05 fm without significantly 
affecting the quality of the mass fits. 

As in the mass formulas Eqs. ([T]) and |2|, one needs two 
parameters bo and co, which control the isoscalar nuclear prop¬ 
erties and thus are needed to reproduce the symmetric nuclear 
binding energy and density. The isoscalar nuclear compressibil¬ 
ity acquires a very reasonable value too, although the saturation 
density is a little low, but not well constrained by the mass fits 
alone. 

Two other parameters a\ and b\ control the symmetry prop¬ 
erties of nuclear matter, and are correlated with a gradient term 
with parameter rj that controls the diffuseness of the nuclear 
surface. These are similar to the parameters «/ and a! 1 in nuclear 

mass formula Eq. (p), but we find that the linear combination 

1 /3 

a\ — b\ p (| is poorly constrained by the masses. This gives 
one an essentially independent control of the “local” symmetry 
energy slope L 2 (not the full L, which is fully determined by the 
neutron equation of state), along with neutron skin thicknesses. 

The addition of quartic isovector coefficients j3 4 permits the 
NEDF to match the neutron matter equation of state without 
significantly affecting the global mass fit. We thus find that the 


• There is a need to consider quartic terms in isospin den¬ 
sity (p„ — p p ) 4 in the NEDF if one aims to describe 
correctly both nuclei and neutron matter within the same 
unified framework, and in particular the neutron star 
crust. The binding energies, charge radii, and neutron 
skin thickness appear to be insensitive to the properties 
of the neutron equation of state, which can essentially be 
fit independently. The position of the neutron dripline ap¬ 
pears to be controlled by the full symmetry energy S and 
its density dependence L, unlike neutron skin thickness, 
which is controlled by the “local" density dependence 
of the symmetry energy L 2 . The properties of nuclear 
matter in stellar environments (when N Z) will there¬ 
fore be controlled by S and L, influencing for example 
the reaction flow in the r-process, the structure of the 
neutron star crust, and the vortex pinning mechanism in 
neutron star crust. 


Entrainment terms av n ■ v p do not appear in any stan¬ 
dard theory of large amplitude collective motion in nu¬ 
clear physics (Baranger and Veneroni| 1978; Brink et «Tj 


1976| Ring and Schuck 2Q04|>, despite being allowed by 


symmetry. They have direct analogues in superfluid mix¬ 
tures - the Andreev-Bashkin effect - and are as natural 
to consider in the presence of mixed proton and neutron 


’He and 4 He superfluids ( 

Andreev and Bashkin 

1975 

Vardanyan and Sedrakyan 

1981; Volovik et al. 1975 ). 


These terms have little influence of the ground state 
binding energies, but a strong effect on excited isovector 
modes and such terms have been discussed in the physics 
of neutron stars and determined to be relevant for the 
description of glitches (Alpar et al. J1984 [Babaev 20041 


|Chamell|2013t|Chamel and Haenself 2006| Gusakov and| 

Haensel||2005l >. 
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To extend the accuracy of this NEDF to reproduce finer details 
of nuclear binding energies (screening of the Coulomb interac¬ 
tion and a correct treatment of the exchange term, isospin and 
charge symmetry breaking effects, Wigner energy, effects due 
to the restoration of various other broken symmetries, etc.) first 
requires an accurate accounting of shell effects; the procedure 
is outlined. Next, the parameters a and £ can be determined by 
fixing the correct value of the Thomas-Reiche-Kuhn sum rule. 
Additional parameters should be introduced into the NEDF 
to model the appropriate physics such as spin-orbit interac¬ 
tions and pairing. The NEDF model is physically intuitive, 
providing a framework for systematically adding physically 
motivated parameters. 

The ability to study dynamical phenomena is a great advan¬ 
tage of the present approach. It is not difficult to develop a 
generalization of the NEDF to finite temperatures and to in¬ 
clude dissipation in the dynamics. In particular it would be 
extremely interesting to establish the general trends in which 
mass and charge asymmetry influences the excitation energies 
and average angular momentum of fragments in induced fis¬ 
sion as a function of the initial energy injected into the nucleus 
and impact parameters in case of reactions. It would be inter¬ 
esting as well to study the sticking probability of two colliding 
heavy-ions as a function of the mass and charge ratio, impact 
parameter and relative velocity. 

A great advantage of the present approach is to provide a 
clear strategy for improving the quality of NEDFs by separating 
the various energy scales. Flere we have identified the bulk 
properties, and shown that they can be properly accounted for 
with a minimal number of parameters. Now we have a clear 
path forward to refine the structure of the NEDF to properly 
account for smaller corrections. In this respect the approach 
outlined here is similar in spirit to an effective field theory. The 
next step is to introduce the spin-orbit coupling, the pairing 
coupling constants (about which there is already an abundance 
of information), and to pinpoint the values of the couplings 
a and |, all of which have a very small effect on the nuclear 
bulk properties. The major impact of these new constants 
will be on the shell-correction energies, and extensive past 
experience indicates that the rms energy should be reduced 


drastically from about 2.6 MeV to close to 0.5 MeV or so (Bohr 


and Mottelson 1998; Brack et al. 1972, 1985:; Moller et al. 

1995 

Moller et al. 2012 

Myers and Swiatecki 

1974 1990 

1991 

1996 

1969 

1966; 

Ring and Sell ucId 2004 

Strutinsky 

1966 

1967 

1968; 

Strutins 

cy and Magnen 1976) or even below. 


accounted for ( 

Aberg 

2002 

Barea et al. 

o 

o 

<N 

Bohigas and 

|Leboeuf| 2002aJb|; Hirsch et al. 

[2005 

2004; 

Molinari and 


IWeidenmullerl [20061 |2004l |Qlofsson efai\ |2006||2008) . 
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VI. SUPPLEMENTAL MATERIAL 

Tables with the results of the numerical fits of various NEDFs 
presented in Table [TlJ along with the Python code used to per¬ 
form these fits and a MATLAB code used to determine the 
ground state and the 3D time-dependent fission process illus¬ 
trated in Fig[l2]can be downloaded from http: //faculty 
Washington.edu/bulgac/NEDF 
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